Surface Granular flows: Two Related Examples 
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Granular surface flows are common in industrial practice and natural systems, however, theoretical 
description of such flows is at present incomplete. Two prototype systems involving surface flow 
are compared: heap formation by pouring at a point and rotating cylinders. Continuum models for 
analysis of these flows are reviewed, and experimental results for quasi-2d systems are presented. 
Experimental results in both systems are well described by continuum models. 
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I. INTRODUCTION 

Granular flows have been the subject of considerable 
recent work [1-5] driven by both technological needs 
(a, u) and the recognition that many aspects of the ba- 
sic physics are poorly understood @|. Surface flows of 
granular materials, that is flows confined to a surface 
layer on a static granular bed, are important in indus- 
trial practice and nature. Industrial examples appear in 
the transportation, processing and storage of materials in 
systems such as rotary kilns, tumbling mixers, and feed- 
ing and discharge of silos. Examples in nature include 
the formation of sand dunes, lava flow, avalanches, and 
transport of sediments in rivers. Although considerable 
progress has been made, theoretical description of surface 
flows is incomplete at present. Several approaches, based 
on different assumptions about the physics of the flows, 
have been proposed [9-19]. A few experimental studies 
are also available [9,19-33]. Most work is focussed on 
two systems: heap flow and rotating cylinder flow shown 
schematically in Fig. 2] 

An important feature of surface granular flows is the 
interchange of particles between the flowing layer and the 
fixed bed. In the case of a rotating cylinder the inter- 
change rate is determined by kinematics since the veloc- 
ity of the fixed bed at the bed-layer interface is known. 
The situation in the case of heap flow is more compli- 
cated. Bouchaud et al. [l^] proposed a phenomenologi- 
cal model (BCRE model) in which the interchange rate 
is determined by the local surface angle. A variation 
of this model proposed by Boutreux et al. [l6| (BRdG 
model) has been broadly validated by continuum models 
|18l [ljj and experiments [l|| , as we show below. Con- 
tinuum models developed previously, for both heaps and 
rotating cylinders, are all based on depth-averaged hy- 
drodynamic equations and differ primarily in the consti- 
tutive equations used. All the models contain parameters 
which must be evaluated from experiments, but in most 
cases, these parameters have not been determined. 



Here we present here a common continuum based 
framework for the analysis of both heap flow and rotating 
cylinder flow. The treatment closely follows that given 
in refs. [jj| and [2i|. Model predictions are compared 
to experimental results and to predictions of previous 
models. The general continuum model is presented first. 
Results for the heap formation and rotating cylinder flow 
are given next followed by conclusions. 



II. GENERAL CONTINUUM MODEL 

Consider a flowing layer on the surface of a granu- 
lar bed assuming the flow is nearly uni-directional in the 
layer and curvature effects are small. The depth averaged 
continuity equation and the x-momentum balance equa- 
tion are simplified using the following assumptions. The 
bulk density in the layer (p) is nearly constant (since the 
dilation of the flowing particles is not too large in the rel- 
atively slow flows being considered). The velocity profile 
in the layer is linear and of the form [31,32] v x = 2u(y/5), 
where u(x,t) is the depth averaged velocity in the layer 
and 5 is the layer thickness. Slow plastic deformation [33] 
is neglected. The shear stress at the interface is taken to 
be 30] 

T xy \y=o = pd 2 f(p)l—^\ - pgS cos [3 tan f3 s (1) 

where d is the particle diameter and tan (3 S is the effective 
coefficient of dynamic friction, with (5 S taken to be the 
static angle of repose. The stress is sensitively dependent 
on the local bulk density and based on recent empirical 
evidence |30j we take f(p) = cS/d with c m 1.5. The 
governing equations then reduce to 
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where T is the flux from the layer into the bed. Further, 
assuming the static friction forces at the heap-layer inter- 
face to be fully mobilized, the Mohr-Coulomb criterion 
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FIG. 1: Schematic view of surface flow systems: (a) Heap flow (b) Rotating cylinder flow, (c) Coordinate system used in the 
analysis. 



yields 



r xy \y=o = - pgd cos (3 tan f3 m , 



(4) 



where tan f3 m is the effective coefficient of static friction. 
Using eq. QJ and the assumptions given above, eq. (0} 
yields 



u = 7J/2, 
with the shear rate given by 

i cos (3 sin(/3 TO - 
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A similar analysis is given by Douady et al. |lq with the 
difference that no stress constitutive equation such as 
eq. (JTJ is used and instead the shear rate in the flowing 
layer is assumed to be constant. 

III. HEAP FORMATION 



where V — g/jcos{3 m . Thus, the continuum model 
yields a source term similar to the BRdG model; the 
scaling of V is also similar to the BRdG model. 

We further simplify the above equations for two dif- 
ferent geometries of heap formation: closed, as shown 
in Fig. ^,, and open in which the end wall (E, Fig. ^Jl) 
is removed. In the open system at steady state, all the 
material entering the system leaves at the far edge of 
the heap and no particles are absorbed or eroded. This 
implies that T = 0, which on substituting into eq. (0 
= Pm = constant. Using these results in eqs. {7J and 
(JSJ) shows that the average velocity (u) and thickness (8) 
of the flowing layer are also constant in open systems. 
The mass flow rate in the system is given by rh = puST, 
where T is the width of the layer. This expression, to- 
gether with eq. (JSJ), gives the following relationship be- 
tween the layer thickness and mass flow rate 



5= [2m/{Tpj)} 1/2 . 



(10) 



Consider a quasi-steady flow (88 /dt, du/dt w 0) and a 
slowly varying interface angle (d[3/dx sw 0) during heap 
formation. The continuity equation (eq. 0) together with 
eq. (O then reduces to 
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and the momentum balance equation (eq. |3J) together 
with eq. Q simplifies to 
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Combining eqs. Q and (JSJ) yields F = gsm([3 m — 
/3)/7COS/3 m , which, for the case when f3 m sw (3, reduces 
to 



T^V(p m -f3), 



(9) 



Experimental results pL9l | based on flow visualization 
studies validate the above predictions, and sample results 
are given below. Fig. |3l shows the variation of the max- 
imum angle of repose with mass flow rate in the system 
for 2 mm steel balls in an open heap system (filled sym- 
bols). The data indicate that m , and thus the coefficient 
of static friction at the heap-layer interface (tan/3 m ), is 
not a constant but increases with the local flow rate. An 
increase in surface angle with flow rate was also reported 
by Lemieux and Durian j29| . Fig. [2) shows the varia- 
tion of the layer thickness (S) with mass flow rate. The 
solid line is a fitted curve of the form 8 oc m 1 ' 2 . This 
indicates agreement with theoretical predictions (eq. I1U|) 
if the product pj is independent of mass flow rate. 

In a closed system (Fig. ^,) , at steady state we must 
have r = constant for the heap to rise uniformly. Inte- 
grating eq. J7J), the layer thickness profile is obtained as 
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FIG. 2: Variation of the (a) maximum angle of repose (/3 m ) and (b) layer thickness (S) with mass flow rate (m) for 2 mm steel 
balls. Filled symbols: open heap system [19]. Open symbols: rotating cylinder system for three cylinder sizes [30]. The solid 
line in (b) is a best fit of the form S oc m 1 ' 2 



5= [81 + 2T(L - x)/j] 
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where 5l is the layer thickness at the end of the layer, 
x = L, and L is the length of the interface (Fig.QJi). The 
rise velocity is related to the mass flow rate by 



T = m/(TLp), 



(12) 



and the interface angle is calculated from eq. JJjJ. 

Experimental results 19] for closed systems show that 
the rise velocity (T) varies nearly linearly with mass flow 
rate in agreement with eq. (|12l) . and the bulk density, 
which is found to be constant, is p — 3.2 g/cm 3 . Fig. [3] 
shows the variation of both interface angle (/?) and layer 
thickness (S) with length along the bed-layer interface 
(L — x) for a fixed mass flow rate. The solid line in 
Fig- EI 3 is a fit of eq. i|ll|) . There is a good match between 
the fitted line and the experimental data, which suggests 
that the shear rate, 7, is constant. Similar results are 
obtained for all flow rates studied. Using experimental 
results for the rise velocity (T) and the interface length 
(L), we obtain 7 = 20 ± 2 s _1 from eq. I|ll|l . where the 
standard deviation indicated is calculated for all 10 flow 
rates studied. Using the value of the bulk density ob- 
tained above, we find from eq. QlOfl that the shear rate 
for the open system is 7 = 22 ± 3 s _1 . The value of 
the shear rate predicted from eq. JBJ is 7 = 20 ± 5 s _1 
for the range of mass flow rates considered. Thus the 
shear rates for the open and closed systems are the same 
within experimental error, and predictions of theory are 
in reasonable agreement with experimental values. 



IV. ROTATING CYLINDER 

The simplest case corresponds to rotating cylinder flow 
for 50% fill fraction. Assuming a nearly flat interface, the 



source term is given by T = ujx. Substituting into the 
continuity equation (eq. |2J) and integrating we obtain 



uS = - (L 2 - x 2 
2 v 



(13) 



Using eq. JSJ the momentum balance equation (eq. |3J) 
simplifies to 



du 3g sin(/3 — f3 s ) 



dx 



cosf3 s 



cdu 2 



(14) 



We consider two different limiting solutions to eqs. 113|) 
and (jT3|) below. 

Firstly, consider the case when shear rate (7) is nearly 
constant. Using eq. JSJ), the flux equation fea. lT^|l gives 
the layer thickness profile as 
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which is symmetric for all rotational speeds (u>). Eq. 1|15|) 
corresponds to the model of Makse 17], in which the 
shear rate is assumed to be a fitting parameter. In the 
present case the shear rate is obtained from eq. © and 
the mean velocity is given by u — jS/2. Substituting 
these results in eq. (|14|l . and using the Mohr-Coulomb 
condition (eq.0J yields eq. © with T = ujx and j3 m ~ /3. 
This allows for calculation of the angle ((3) along the in- 
terface. Thus the assumption of a constant shear rate is 
consistent with the model, and gives a complete descrip- 
tion of the flow. However, it is not apparent from the 
analysis, under what conditions the solution is valid. 

Consider next the case when the shear rate is not con- 
stant along the layer, but when the acceleration (du/dx) 
is small. Eliminating S using eq. (|13f) . the scaled momen- 
tum balance becomes 
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FIG. 3: Variation of the (a) surface angle (/3) and (b) layer thickness (/3) with distance from the edge of the heap (L — x) for 
flow of 2 mm steel balls in a closed system. Symbols are experimental data and error bars indicate the standard deviation over 
six measurements. Solid line in (a) is a fit of eq. © and in (b) is the prediction of eq. lilt . 



where u — u/luL, £ = x/L and the dimensionless param- 
eters are the Froude number, Fr — u> 2 L/g, and the size 
ratio, s = d/L. The first term on the right hand side 
of eq. i|16fl is the net driving force, that is the gravita- 
tional force less the frictional resistance to flow, and is 
independent of the flow velocity (u). The second term 
is the 'viscous' resistance due to collisional stresses, and 
the third term arises as a result of in-flow and out-flow 
of particles from the layer. Both these terms depend on 
the flow velocity. Typical experimental Froude numbers 
for experiments in rotating cylinders are in the range 
O(10~ 3 ) to 0(1O~ 2 ). In these cases the driving force 
term (0(1/ Fr)) is much larger than the acceleration term 
(0(£/VsFr) based on eq. 0), particularly near the mid- 
point of the layer (£ = 0) . The collisional stress term is 
of the same magnitude as the net driving force term since 
the flow velocity increases to balance the two. Thus for 
^^Fr/s <C 1 the acceleration term may be neglected. 

For negligible acceleration (du/d£ ~ 0), the scaled 
mean flow velocity is obtained from eq. (|1(JH as 
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where A = sin(/3 — /3 S )/ cos (3 S 
layer thickness profile is 



Using eq. ifH?)) . the scaled 
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(18) 



where 6 = S/L. The above solution is valid only if A > 0, 
that is if j3 > (3 S . For (3 < (3 S , we have u = 8 = 0, thus 
there is no steady flow possible if the interface angle is 
less the static angle of repose. This is consistent with 
the definition of the static angle of repose. Note that the 
layer profile is not symmetric about £ = 0, and for any 
£ > we have 6(—£) > S(Q, that is, the upper part of 
the layer (£ < 0) is thicker than the lower part. The 



source of the asymmetry is the in-flow/out-flow term in 
the momentum balance (third term on the right hand 
side of eq. Hfl|l . In the upper part of the layer (£ < 0) 
the flow is retarded by material entering the layer from 
the bed (T < 0) and the reverse is true in the lower part 
of the layer. Thus, the layer is thicker in upper part be- 
cause of the lower velocity relative to the lower part of 
the layer (£ > 0), resulting in a skewed profile. Further, 
eq. (|18fl indicates that the profile becomes more skewed 
with increasing Froude number (Fr) and decreasing size 
ratio (s). In the limit, Fr/s <C 1, the scaled layer thick- 
ness profile becomes 8 = (csFr/A) 1 / 4 ^ — £ 2 ) 1 ^ 2 , which 
is identical to the result obtained assuming a constant 
shear (eq. Il5|l when eq. © is used to calculate the shear 
rate. This implies that a profile symmetric about the 
layer midpoint (£ = 0) is obtained at very low Froude 
numbers and relatively high size ratios, and in this limit 
the shear rate is nearly constant. 

The interface angle profile (/?(£)) is obtained from 
eq. |JSJ, using T — lux and eq. iJBJ as 



(17) m = Pr, 



Fr cos /3„ 



3cs 



£ + (£ 2 + 9csA/Fr) 



1/2 



1/2 



(19) 

In simplifying the preceding equation we assume (3 m *w 
P w j3 8 . Eq. (|19fl indicates that the interface angle de- 
creases monotonically with distance along the interface 
and at £ = 0, /3 = (3 m . Thus in the rotating cylinder 
flow the maximum angle of repose can be experimen- 
tally obtained by measuring the interface angle at the 
midpoint of the layer. For (Fr/s)£, sufficiently large and 
£ > 0, we get /3 < 0, that is, for small size ratios and 
large Froude numbers the layer profile may turn up at 
the end. Conversely, when (Fr/s) <C 1, eq. lfH?)l yields 
(3 w f3 m , and the interface profile is nearly flat. Ne- 
glecting terms 0(£Fr/s), which is consistent with the 
approximation in the momentum balance equation, we 
get A = (I3 rn - I3 S )/ cos/3 s . 
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FIG. 4: Variation of the layer thickness at the midpoint (5(0)) with (sFr/A) ' for (a) steel balls, (b) glass beads (c) and sand 
particles in cylinders of different sizes and at different rotational speeds. Symbols are experimental data for different sized 
particles: o d = f mm, Ad = 2 mm, O d — 4 mm, V d — 0.4 mm and X d = 0.8 mm. The solid line is a fit of eq. (1181 and the 
values of the parameter c ' are indicated. 
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FIG. 5: Layer thickness profiles for (a) 2 mm steel balls and (b) 0.8 mm sand. Symbols denote experimental data for three 
different froude numbers (Fr): o Fr — 2 x 10 -3 , A Fr — 22 x 10~ 3 , O Fr = 64 x 10~ 3 . Solid lines are predictions of eq. itTSt 
and dashed lines are the predictions of the model of Khakhar et al. [13]. The error bars give the standard deviation over 6 
measurements and the bar indicates the scaled diameter of a particle (s = d/R). 



Consider next a comparison of the theoretical results to 
experimental data. A few key numbers are reported, as 
they convey a sense of qualitative agreement. However, 
for full details the reader is referred to J3(j. The model 
parameters required are (3 S , (3 m and c. Data of Orpe and 
Khakhar 30] for the first two parameters are shown in 
Fig. [2^ (open symbols) for 2 mm steel balls in rotating 
cylinders of 3 sizes and for different rotational speeds of 
the cylinders. The data correlates reasonably well with 
the mass flow rate at the midpoint of the layer calculated 
from rh = pu>L 2 T/2, where T is the cylinder length and 
the same density as in the heap experiments (p = 3.2 
g/cm 3 ) is used. Data spanning nearly two decades of 
flow rate fall on a single curve, although with some scat- 
ter. The maximum angle of repose increases with mass 
flow rate, and the measured values are similar to those 
from heap experiments which are also shown in the same 



figure. The static angle of repose is the angle at m = 0. 

Orpe and Khakhar [3(| had obtained c ~ 1.5 by fit- 
ting the theory of Khakhar et al. [Tij to experimental 
layer thickness profiles. We obtain a new estimate of 
the parameter based on the layer thickness at the mid- 
point of the layer (£ = 0), which, from eq. I|18[) . is 
(5(0) = (csFr/A) 1 / 4 . Fig. Q] shows experimental data for 
6(0) versus (sFr/A) 1 / 4 for experimental data for 90 ex- 
periments comprising steel balls, glass beads and sand 
of different sizes in cylinders of different sizes and for 
different rotational speeds. The data falls on a straight 
line for each material (although with some scatter) and 
a least squares fit gives c = 1.9 for steel balls, c = 1.6 
for glass beads and c = 1.4 for sand. Since the model 
is essentially exact at £ = 0, the good fit implies that 
the proposed constitutive equation for stress is reason- 
able, and the shear rate in the layer is well-described by 



eq. © at £ = 0. 

Predictions of the model for the layer thickness pro- 
file are compared to experimental data in Fig. for sand 
particles and steel balls for different Froude numbers in 
a cylinder of radius 16 cm, using the value of c obtained 
above and experimental values for (3 m and j3 s . The agree- 
ment is good except at the highest Fr and low s, and the 
all the qualitative features of the data are reproduced. At 
low Fr and relatively high s studied, the profile is nearly 
symmetric (steel balls at the lowest Fr), and the profiles 
become more skewed with increasing Fr and decreasing 
s. The deviation at the high Froude numbers and low size 
ratio are due to neglect of the acceleration term. Similar 
agreement is obtained for the other cases studied as well. 
The predictions of the model of Khakhar et al. [lj| are 
shown in the figure as dashed lines. These nearly coincide 
with the results from the present model, except for the 
highest Fr for sand, indicating that the approximations 
made are reasonable for the parameter values of interest. 
It is remarkable that such a simple theory is able to de- 
scribe the behavior of the system over such a wide range 
of parameters: materials include steel balls, glass beads 
and sand; varying shapes with steel balls being spherical, 
glass beads, nearly spherical and sand being irregularly 
shaped; size ratios in the range s £ (0.005, 0.05) and 
Froude numbers in the range Fr £ (2 x 10~ 3 , 64 x 10~ 3 ). 
Model predictions of the interface angle profile are in rea- 
sonable agreement with experiments 30]. 



form thickness with a uniform flow velocity is obtained, 
whereas in the closed heap system 8 2 ex x. The inter- 
face angle is constant and equal to the maximum an- 
gle of repose in the open system, whereas it decreases 
with distance from the pouring point in the closed sys- 
tem. Results for the rotating cylinder are obtained for 
the case when the acceleration of particles in the layer 
is small ($,\/Fr/s <C 1). The layer profile is found to 
be asymmetric about the midpoint of the layer (£ = 0) 
with the upper part of the layer (£ < 0) being thicker. 
The skewness increases with increasing Froude numbers 
and decreasing size ratios. The scaled shear rate (j/ui) 
decreases with increasing Froude number and size ratio. 
The layer interface angle decreases with distance in the 
flow direction. For high t;Fr/s and £ > the layer turns 
up, whereas when £Fr/s is small a nearly flat interface 
is obtained. 

Quasi-2d experiments carried out for open and closed 
heaps and rotating cylinders of different sizes, by and 
large, validate the predictions of the theory. The three 
material parameters of the model (j3 s , (3 m and c) are 
all obtained from relatively simple measurements. The 
model equation can thus be applied to more complex 
geometries. Deviations of the model from experimental 
data appear in the interface angle profile in the rotating 
cylinder flow. This is most likely due to end wall effects 
which are discussed in 301. 



V. CONCLUSIONS 

A theoretical framework serves to unify the behaviour 
of surface flows for two prototypical systems: heap flow 
and rotating cylinder flow. The model is based on a 
stress constitutive equation and failure criterion which 
contain three material parameters: /3 S , f3 m and c. Ana- 
lytical results for both systems give a complete descrip- 
tion of the systems in terms of the layer thickness profiles 
(8(x)), average velocity of flow (u(x)) and the interface 
angle profile (j3(x)). In open heap systems a layer of uni- 
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